Objectives

Download NHD layers for Klamath Basin

Load NHD layers from geodatabase

layers <- st_layers(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'))

nhd_fcode <- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'NHDFCode') |> 
  janitor::clean_names() 
## Reading layer `NHDFCode' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
nhd_flowline_vaa<- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'NHDPlusFlowlineVAA') |> 
  janitor::clean_names() 
## Reading layer `NHDPlusFlowlineVAA' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
nhd_flowline_raw <- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'NHDFlowline') |> 
  janitor::clean_names() 
## Reading layer `NHDFlowline' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 341762 features and 18 fields
## Geometry type: MULTILINESTRING
## Dimension:     XYZM
## Bounding box:  xmin: -124.4217 ymin: 38.29819 xmax: -120.5703 ymax: 43.33625
## z_range:       zmin: 0 zmax: 1612
## m_range:       mmin: 0 mmax: 100
## Geodetic CRS:  NAD83
huc_12 <- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'WBDHU12') |> 
  janitor::clean_names() |> 
  st_transform(crs = "+proj=longlat +datum=WGS84") |> 
  filter(substr(huc12, 1, 6) == "180102")
## Reading layer `WBDHU12' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 705 features and 21 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.5351 ymin: 38.25486 xmax: -120.5637 ymax: 43.34273
## Geodetic CRS:  NAD83
huc_10 <- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'WBDHU10') |> 
  janitor::clean_names() |> 
  st_transform(crs = "+proj=longlat +datum=WGS84") |> 
  filter(substr(huc10, 1, 6) == "180102")
## Reading layer `WBDHU10' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 147 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.5351 ymin: 38.25486 xmax: -120.5637 ymax: 43.34273
## Geodetic CRS:  NAD83
huc_8 <- st_read(here::here('data-raw', 'shapefiles', 'NHDPLUS_H_1801_HU4_GDB.gdb'), layer = 'WBDHU8') |> 
  janitor::clean_names() |> 
  st_transform(crs = "+proj=longlat +datum=WGS84") |> 
  filter(substr(huc8, 1, 6) == "180102")
## Reading layer `WBDHU8' from data source 
##   `/Users/maddeerubenson/Documents/git/rivermmile/data-raw/shapefiles/NHDPLUS_H_1801_HU4_GDB.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 22 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.5351 ymin: 38.25486 xmax: -120.5637 ymax: 43.34273
## Geodetic CRS:  NAD83
fcode_descr <- c("Stream/River", "Artificial Path",
                 "Stream/River: Hydrographic Category = Intermittent",
                 "Stream/River: Hydrographic Category = Perennial", "Stream/River: Hydrographic Category = Ephemeral")  

nhd_flowline <- nhd_flowline_raw |> 
  mutate(rc_huc_6 = substr(reach_code, 1, 6),
         rc_huc_8 = substr(reach_code, 1, 8),
         rc_huc_10 = substr(reach_code, 1, 10),
         rc_huc_12 = substr(reach_code, 1, 12)) |> 
  filter(rc_huc_6 == "180102") |> # filter to klamath huc6
  filter(!is.na(gnis_name)) |> # remove NA streams 
  left_join(nhd_fcode) |> 
  filter(description %in% fcode_descr) |> 
  left_join(nhd_flowline_vaa) |> 
  glimpse()
## Rows: 31,796
## Columns: 73
## $ permanent_identifier         <chr> "137190692", "147512448", "122413487", "8…
## $ f_date                       <dttm> 2012-02-26 18:16:31, 2012-02-26 18:23:10…
## $ resolution                   <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
## $ gnis_id                      <chr> "01141109", "01156651", "01657346", "0112…
## $ gnis_name                    <chr> "Doe Creek", "Long Creek", "Scaath Creek"…
## $ length_km                    <dbl> 3.35500000, 5.33000000, 0.08269511, 4.451…
## $ reach_code                   <chr> "18010201000166", "18010202000293", "1801…
## $ flow_dir                     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ wb_area_permanent_identifier <chr> NA, NA, "122415943", "84069699", NA, NA, …
## $ f_type                       <int> 460, 460, 558, 558, 460, 460, 460, 460, 4…
## $ f_code                       <int> 46003, 46006, 55800, 55800, 46003, 46003,…
## $ main_path                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ in_network                   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ visibility_filter            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ shape_length                 <dbl> 0.0374853923, 0.0606200004, 0.0007596266,…
## $ nhd_plus_id                  <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ vpuid                        <chr> "1801", "1801", "1801", "1801", "1801", "…
## $ enabled                      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ rc_huc_6                     <chr> "180102", "180102", "180102", "180102", "…
## $ rc_huc_8                     <chr> "18010201", "18010202", "18010209", "1801…
## $ rc_huc_10                    <chr> "1801020100", "1801020200", "1801020900",…
## $ rc_huc_12                    <chr> "180102010001", "180102020002", "18010209…
## $ description                  <chr> "Stream/River: Hydrographic Category = In…
## $ canal_ditch_type             <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ construction_material        <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ hydrographic_category        <chr> "Intermittent", "Perennial", " ", " ", "I…
## $ inundation_control_status    <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ operational_status           <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ pipeline_type                <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ positional_accuracy          <chr> NA, NA, " ", " ", NA, NA, NA, NA, NA, NA,…
## $ relationship_to_surface      <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ reservoir_type               <chr> " ", " ", " ", " ", " ", " ", " ", " ", "…
## $ stage                        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ stream_leve                  <int> 4, 4, 5, 4, 4, 4, 7, 7, 5, 5, 5, 7, 4, 6,…
## $ stream_orde                  <int> 2, 4, 2, 6, 2, 2, 3, 3, 2, 3, 2, 3, 3, 1,…
## $ stream_calc                  <int> 2, 4, 2, 6, 2, 2, 3, 3, 2, 3, 2, 3, 3, 1,…
## $ from_node                    <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ to_node                      <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ hydro_seq                    <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ level_path_i                 <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ path_length                  <dbl> 0.00000, 16.63400, 13.92796, 0.01400, 17.…
## $ terminal_pa                  <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ arbolate_su                  <dbl> 8.332000, 73.236400, 3.444947, 10043.2644…
## $ divergence                   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ start_flag                   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ terminal_fl                  <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ up_level_pat                 <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ up_hydro_seq                 <dbl> 5.00004e+13, 5.00004e+13, 5.00004e+13, 5.…
## $ dn_level                     <int> -1, 4, 4, 4, 4, 4, 7, 6, 5, 5, 4, 7, 4, 5…
## $ dn_level_pat                 <dbl> 0.00000e+00, 5.00004e+13, 5.00004e+13, 5.…
## $ dn_hydro_seq                 <dbl> 0.00000e+00, 5.00004e+13, 5.00004e+13, 5.…
## $ dn_minor_hyd                 <dbl> 0.00000e+00, 5.00004e+13, 0.00000e+00, 5.…
## $ dn_drain_cou                 <int> 0, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ from_meas                    <dbl> 0.00000, 0.00000, 0.00000, 0.00000, 0.000…
## $ to_meas                      <dbl> 59.11259, 100.00000, 5.04551, 40.59847, 1…
## $ rtn_div                      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ thinner                      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ vpu_in                       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ vpu_out                      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ area_sq_km                   <dbl> 8.67549996, 9.59290002, 0.00760000, 3.879…
## $ tot_da_sq_km                 <dbl> 12.8034, 87.4999, 1.6544, 10137.9325, 34.…
## $ div_da_sq_km                 <dbl> 12.8034, 87.4999, 1.6544, 2510.5202, 34.3…
## $ max_elev_raw                 <dbl> -9998, -9998, -9998, -9998, -9998, -9998,…
## $ min_elev_raw                 <dbl> 140977, 155705, 655, 124234, 149933, 1489…
## $ max_elev_smo                 <dbl> 147172, 165024, 800, 124607, 153047, 1644…
## $ min_elev_smo                 <dbl> 140977, 155705, 655, 124252, 149977, 1489…
## $ slope                        <dbl> 0.01846498, 0.01748405, 0.01753429, 0.000…
## $ slope_len_km                 <dbl> 3.35500000, 5.33000000, 0.08269511, 4.451…
## $ elev_fixed                   <int> 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1,…
## $ hw_type                      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ hw_node_sq_km                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ status_flag                  <chr> "A", "A", "A", "A", "A", "A", "A", "A", "…
## $ Shape                        <MULTILINESTRING [°]> MULTILINESTRING ZM ((-121…
ggplot() +
  geom_sf(data = hucs, aes(color = 'darkred')) +
  geom_sf(data = nhd_flowline) +
  theme_minimal() +
  theme(legend.position = "none")

Filter NHD layer to mainstem rivers/tribs

# filtering to mainstem and tribs of interest 

rivers <- c('Klamath River', 'Trinity River', 'Sprague River',
            'Williamson River', 'Wood River', 'Lost River', 'Shasta River', 
            'Scott River', 'Salmon River', 'Blue Creek', 'Bogus Creek', 'Clear Creek',
            'Indian Creek', 'Link River', 'Jenny Creek')

nhd_flowline_filtered <- nhd_flowline |> 
  filter(gnis_name %in% rivers) |> 
  st_transform(crs = "+proj=longlat +datum=WGS84") |> 
  st_zm()

View filtered NHD layers

leaflet(data = nhd_flowline_filtered) |> 
  addTiles() |> 
  addPolylines(popup = ~paste0("name: ", gnis_name, "<br>",
                               "id: ", gnis_id))

Clean NHD layer

I decided to process the NHD file in QGIS - the project is saved within the data-raw/shapefiles folder.

Processing Notes:

# Read in processed 

processed_nhd <- sf::read_sf(here::here('data-raw', 'shapefiles', 'cleaned_nhd_layer.shp'))

leaflet(processed_nhd) |> 
  addTiles() |> 
 addPolylines(popup = ~paste0("name: ", processed_nhd$gnis_nm, "<br>",
                               "id: ", processed_nhd$gnis_id))

Collapse NHD segments into single lines

klamath_rivers <- processed_nhd |> 
  group_by(gnis_nm) |> 
  summarise(geometry = st_union(geometry)) |>  
  st_cast("MULTILINESTRING") 
 # st_transform(crs = 32610) 

leaflet(klamath_rivers) |> 
  addTiles() |> 
  addPolylines(popup = ~paste0(gnis_nm))

Save cleaned rivers file

usethis::use_data(klamath_rivers, overwrite = TRUE)
## ✔ Setting active project to "/Users/maddeerubenson/Documents/git/rivermmile".
## ✔ Saving "klamath_rivers" to "data/klamath_rivers.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
knitr::knit_exit()